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Codes involving one and two spatial dimensions and three velocity 
dimensions are used to model the earth's magnetotail. It is shown that the 
magnetotail can become Inflated as a consequence of low energy plasma 
convection toward the neutral plane. The two-dimensional code shows the 
development of small-scale turbulence. However, the growth of the turbulence 
appears limited, and does not lead to substantial dissipation. A one- 
dimensional model was used to simulate a particle population consisting of 
plasma sheet ions and cold lobe plasma. An applied convection electric field 
causes plasma sheet thinning and later expansion, consistent with observed 
substorm morphology. During the course of a simulated substorm the cold 
particles are heated and merge with the plasma sheet population. The magnetic 
field energy shows a slight increase, followed by a more rapid decrease as all 
of the cold particles are swept into the magnetotail. The code exhibits a 
conversion of both magnetic field energy and of energy supplied by the 
convection electric field into particle energy. The simulations suggest that 
much of the magnetotail substorm morphology may be a simple consequence of an 
increase, followed by a decrease, in the convection electric field, without 
the requirement of any magnetospherie size scale plasma instability or other 
disruptive processes. It is also concluded that the presence of the 
convection electric field and a continuing replenishment of low energy 
particles in the magnetotail lobes are both necessary for maintenance of the 
magnetotai 1 . 


1. Introduction. 


Ever since the discovery of the neutral and plasma sheets ln the earth's 
magnetotall considerable attention has been focused upon the relationships 
between the morphology of these features and the dynamics of auroral 
substorms. For a comprehensive review of this subject, see Akasofu ( 1979 ). 
It Is generally accepted that (Akasofu, 1979 ; Fairfield et al., 1981 ) that the 
onset of a substorm, as determined from ground-based Indicators, coincides 
with the thinning of the plasma sheet followed by a subsequent expansion. The 
substorm morphology of the magnetic field In the lobe does not seem to show as 
constant a pattern of Increase or decrease (Akasofu, 1979 ). A southward 
turning of the field In the magnetotall Is sometimes seen, which has been 
Interpreted as a consequence of the formation of an x-type neutral line 
earthward of the point of observation (Hone?, 1979 ). Others (Akasofu, 1979, 
and references contained therein) have Interpreted the southward turning as 
simply the tilting of the magnetic field. 

One Intriguing feature of the earth's magnetotall 1s that 1t contains a 
large amount of free magnetic energy that could be converted to particle 
energy If the plasma sheet dynamics were to permit a relaxation to a more 
"dipolar" configuration. Schindler ( 1974 ) has suggested that the plasma sheet 
thinning during the growth phase of a substorm leads to an ion tearing mode, 
instability which results in a release of magnetic energy. Birn ( 1980 ) has 
developed a numerical fluid inodel for the magnetotall. This model Indicates 
stability if there is no viscous or resistive dissipation In the plasma; but a 
conversion of magnetic or kinetic energy does occur when the model contains 
dissipation. 

Akasofu (1980a, b), on the other hand has Identified a solar wind 
parameter, e, Involving the solar wind velocity and interplanetary magnetic 
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field, with the dimensions of energy flux, that shows a remarkable correlation 
with the magnetospheric substorm AE Index. The fact that the c and AE 
parameters have such a similar variation suggests that the substorm Is 
externally driven, rather than the results of an Instability Internal to the 
magnetosphere (Swift, 1979). 

For this reason, this paper adopts a very different approach to 
magnetotaH dynamics by allowing energy Input Into the magnetotall by a 
convection electric field, and Investigating processes that may lead to the 
formation of the highly stretched tall. The model also Includes processes for 
energy exchange between the tall field and particles, so that processes 
investigated by Schindler (1974) and Birn (1980) can occur In the model to be 
described below. However, since a particle model will be used, there are no 
dissipation processes other than that allowed by the equations of motion of a 
charged particle in the electromagnetic field. Moreover, the model does not 
describe equilibrium processes or departures from equilibrium so the results 
are not easily compared to results of stability analysis. 

The approach taken follows along some of the Ideas expressed by Cowley’ 
(1980). He suggests that very low energy plasma, of solar wind or Ionospheric 
origin finds Its way Into the lobes of the earth's magnetotail, whence it is 
swept into the neutral sheet by the convection electric field. In the 
nonadiabatic encounter with the neutral sheet, the plasma Is accelerated 
earthward and heated (Swift, 1977), thus converting cold lobe plasma into a 
population with the characteristics of the plasma sheet. These ideas are 
given further support by the recent observations with the Hawkeye satellite. 
Anderson and Eastman (private communication) have made measurements of the 
plasma line cutoff in the lobes of the iregnetotail which Indicate plasma 
densities In the range of 0.1 to 1.2 cm‘^. Since the presence of these 
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particles is not indicated by the particle detectors, their energy must be 
below detection threshold, hence they are likely of ionospheric origin. * 

In summary, the calculations will be directed toward answering the 
question of whether the main characteristics of the observed tail morphology 
can be understood in terms of a process of polar wind, or magnetosheath ions 
being swept from the lobes into the neutral sheet region by a convection 
electric field. To this end, we shall consider both one and two-dimensional 
models in a plane continuing the magnetic field, which is assumed to be the 
magnet ospheric noon-midnight meridian plane. The two-dimensional code can be 
used to investigate possible development of turbulence and 0- and X-type 
neutral lines (Stern, 1979; Vasyliunas, 1980). The one-dimensional code can 
be used, because of the much greater economies of operation, to observe the 
system over a much longer timi span in relation to the particle gyroperiod. 
It turns out that turbulence does develop in the two-dimensional simulation, 
but that the turbulence does not seem to affect the large-scale morphology. 
As a result, the one-dimensional model appears to offer a good representation 
of the total dynamics. 

In the next section, the model will be described, along with a brief 
description of th& numerical procedures. The following section will present 
the results of the calculations, and the final section will present a 
comparison between the model results and the observed substorm morphology of 
the earth's magnetotail. 

2. Calculation Procedure. 

2.1 Formulation of the Model. 

The basic code used is a nonradiative electromagnetic code similar to 
that described by Nielson and Lewis (1976) in which we neglect all 
electrostatic Interactions among particles and assume uniformity in the y- 








direction. Electrostatic Interactions are neglected because we also assume 
the existence of a highly mobile population of cold electrons which 1s capable 
of shorting out potentials developed along field lines.. Possible 
electrostatic effects will be described In the concluding section. We further 
assume that the magnetic field Is confined to the x-z planes. This and the 
assumption of y-symmetry Implies that the only allowed currents are In the y- 
dlrection. This also means that currents carried by E x 8 and field-aligned 
motion of Ions and electrons cancel and that the x and z-components of the 
current due to gyratlonal motion and polarization drifts of the Ions average 
to zero. The assumption of y-symmetry also precludes Instabilities excited 
by cross-tall current flow (Huba et al., 1978; Tanaka et al., 1981) which may 
lead to anomalous resistivity and electron heating. 

These assumptions make It possible. to describe the entire system In terms 
of the y-component of the vector potential of which the only source is the y- 
cotiponent of the ion motion. We therefore consider the following system of 
equations for the ion motion. 
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where Ey is the electric field, which includes both a convection electic field 
and the solenoidal field due to the time rate of change of 5. The magnetic 
field ^ is expressed in terms of the vector potential 
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where Ay is the y-component of the vector potential, which satisfies the 
equation 

V^Ay{r) - - 4weE S(r - r,^) Vy|^ (3 

where the sum is over all particles present in the system, and is the 
instantaneous position of the particle, and Vy|^ its velocity along y. S(r) is 
a normalized density function which differs from a simple 6-function because 
in a numerical simulation the charge must be assigned to neighboring grid 
points, is found by integrating 


with V calculated from (1). In (1) the field terms are actually evaluated 
according to 


5 (rj^) (r) S (r - r^) dxdz (S 

At this stage it is useful to introduce dimensionless units in order that 
the system can be scaled to run with economically practical particle numbers, 
and distance and time scales. To this end, we introduce dimensionless 
position and time coordinates. 

t = r/p (6a 
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T » Ut 
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where p and m are a nominal gyroradius and gyrof requency . The dimensionless 
velocity becomes 

V * - V/(lDp) 


The equations of motion (1) become 


dv 
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(c + V X 3)/a. 
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where ^ = e&/(mc), and e 1s the electric field acceleration e*ef/(mti)p) . A 
dimensionless vector potential is introduced 
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such that 
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with these substitutions, the equation for Ky is given by 
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where a is the coupling constant 
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and 


y - p^E s (I - f^) y (12 

k 

The problem is scaled for computation by selecting the number of particles and 
the size of the domain, and then selecting an appropriate value for the 
coupling constant. We can also, without loss of generality, set u • 1, and 
for the sake of convenience make the variable replacements 
I X, V V, Xy + Ay. It will also be useful to split the electric field 

acceleration Into the solenoidal, ■ 5Ay/8t, and the convective term 

which represents the dawn-to-dusk convection electric field. 

The particle and field equations In the new dimensionless , simulation 
units, now read 


dv 

IT • ""z^ 

(13a 

dVj 

3T • ’“/y 

(13b 

dP 

(13c 


where the field variables at the position of the kth particle are computed 
according to the prescription given in (5). Py is the canonical momentisn and 
is used to calculate Vy. 

Vyk = Pyk ’ ^ ' 
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Use of the canonical momentun makes It possible to avoid the explicit 
appearance of the destabilizing aA^/dt term In the equation for Vy (Nielson 
and Lewis, 1976). The equation for the vector potential, can be written 
In the form 


7^Ay(r) - -a(P(r) + n(r)Ay(r)) (15 

where 

P(r) E Py,^ S (r-r^) 
and 

n(r) • E S (r-ri.) 
k ^ 

The particle gyrofrequency Is calculated from 

6A 

(®x» “z^ * 

In addition to the field generated by the local particle currents, we 
will wish to include a constant, uniform magnetic field in the z-directlon. 
The earth's dipole field, and possibly the Interplanetary magnetic field, 
could* contribute a z-component, but not necessarily uniform field. However 

making provision for a nonuniform, asymptotic z-component, would greatly 
complicate the calculation. Hence, Ay in (15) will Include an additive term. 
Ay Ay + Q^x , which gives rise to a uniform magnetic field in the z- 
directlon of intensity 
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(16a 

(16b 




Before describing the numerical procedure. It will be useful to make note 
of important relations expressed in the simulation units of (13)-(17). 
Multiplying the equations (13a, b) by and respectively, and the equation 

a e-v Q -*■ V fl 
X z z X 

by Vy and summing over all the particles, the energy conservation law is 
obtained 



^ C Z J dxdz] - / e fy X S • ds (18 

where the integral over h is an integral over a boundary surface. From this 
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expression we identify 0 /2a as the magnetic energy density and 

ety X Q/a SIS the Poynting flux. The E x B drift is Ug.»Cj.ty x S/p^. Another 

important parameter is the Alfvdn velocity, V^, which given by the expression, 

* a^/(on) (19 

a 

where n is computed from (16b). 

2.2 Two-dimensional Numerical Procedures. 

The integration of the particle equations of motion (13a, b) and the non- 
dimensional version (4) is accomplished by the standard leap-frog scheme, with 
the X- and z-components of the velocities and positions being advanced a half 
time-step apart. The value of Py in (13c) is advanced following the 

calculation of Ay; thus Vy, which is calculated from (14), is at the same 
time-step of the position coordinates, but a half time-step removed from the 


(Vx» Vy) variables. Therefore, the particle pushing algoritivn is second order 
accurate. The integration of (13c) is only first order accurate, but since Py 
changes nearly linearly with time this introduces little error. 

An iteration scheme is used to solve (14) is similar to that used by Lin 
(1978), namely 

7^Ay*'*‘^ - on Ay^^^ ■ -oP + o(n-n) (20 

where the superscript l refers to the iteration number. The A«0 iteration is 
the last iteration of the previous time-step, n is a constant parameter. If 
fi is too small, the iteration is unstable, and if n is too large, the 
convergence will be slow. This slow convergence is due to the fact that the 
Green's function for the operator en the right hand side is 
G * exp(-Ry4n)/R, where R»|r-r |. Thus, if there is a change in the current 
at r , the effect of this change wi 1 1 . propagate only a distance (on)“ each 
iteration. Large values of n require either a large number of iterations 
between time-steps, or very short time-steps in comparison to the time scale 
of the simulation, in order that changes in one part of the domain might 
propagate across the entire domain. Nielson and Lewis (1976) suggest 
" “V2(n^ax ^ "min)- 

The equation for A *+i is a constant -coefficient, Helmholtz equation of 
the type 

(7^-k^) A(r) » -p(r) (21 

In solving this equetion, we can neglect the Q^x term of the constant 
background field, because in the numerical procedure the term anQ^x is 
incorporated into the aP term on the right-hand side of (20). This makes it 
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possible to solve (21) subject to periodic boundary conditions in x for a 
domain of length L. 

Equation (21) is solved using Fourier transform methods. The details of 
the method which allow the Imposition of boundary conditions at z « 0 and z ■ 
L is given in the appendix. Essentially, the boundary conditions are that the 
X'Component of the magnetic field, averaged over x at z > 0 and L be equal and 
opposite and that the similarly averaged boundary values of the vector 
potential be constant. 

The source terms in the equation of /^, shown in {16) are confuted by a 
standard four-point area-weighting scheme in which the contribution of the 
current is distributed at the four surrounding grid points. Momentum 
conservation therefore requires that the interpolation of the fields as 
indicated in (5) be the usual bilinear four-point Lagrangian scheme. 
Therefore, each particle contribution to the current is represented as values 
at the four grid points defining the cell which contains the particle. No 
other smoothing of the particle density such as the Gaussian weighting, is 
used, because this would introduce inconsistencies in the boundary conditions. 

The particle boundary conditions are assumed periodic in x. That is if 
X > L, the replacement x x-t is used, and Vj^ and v^ are unchanged. A 

similar replacement is used if x < 0. However, because of the background 

magnetic field, Q^, the value of Py must be changed by an amount iQ^l when a 
particle crosses the x * 0, L boundaries. It is this change in Py on crossing 

the X = 0, L boundaries that would give rise to problems if a Gaussian charge 

weighting scheme were used. If a particle were located at x * L - 6, part of 
the charge would appear at the x * 0 grid points, with a totally inappropriate 
value of Py attached, which would give rise to spurious boundary currents. 
When the particle crosses the z * G. L boundaries, which is an unusual event 
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because of the inward convection, it is reinserted in the doinain at the same 
boundary, on the same guiding center, but with the velocity reinitialized to 
the E X B drift. 

2.3 One-Dimensional Numerical Procedures. 

In the one-dimensional problem all field variables depend only on the z- 
coordinate, except for the vector potential which is written 

Ay(x, z) ■ a(z) + Q^x (22 

where is constant. The field equation for a is therefore 

. .a ^ . a(z) * S (r - r^) (23 

Upon integrating (23) over x frorn 0 to L and dividing by L, we obtain 

& - ® n(z)a - - -S r (p , 9 x.)S(z - z.) (24 

dz^ L L 1 ^ yK z K K 

where 

n(z) - I S(z-zJ (25 

k ^ 

since the only x-dependenee in (24) occurs with Py|^, we can define a new 
variable 




“2*k 


and (13c) is replaced by 


(26 
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dt 




where is advanced by (13a), and It is no longer necessary to keep track of 
for each particle. The field equation is now given by 


where 


2 

^ - (x^n(z)a » -o^P(r) 
dz 

P - I ?y„S{z - z,,) 


(28 

(29 


and is the one-dimensional charge coupling constant. V and n are 

calculated using a linear weighting scheme and the field values at the 
particles are calculated from grid points by linear interpolation. 

No iteration is necessary in solving (28) because it can be solved 

efficiently for variable n(z) by using the finite difference approximation to 
the derivative, and using the tri -diagonal algorithm (Roache, 1976) for 
calculating a(z) at the grid points. 

In the problem at hand, we can assume reflection symmetry about the z = 0 
plane. It is therefore sufficient to consider only the domain from the center 

of the neutral sheet outward into the lobe of the magnetotail. Since the x- 

component of the magnetic field vanishes, one boundary condition is that 
aa/bz = 0 at z a 0, the assumed position of the neutral sheet. At the other 
end of the domain, we wish the value of a to be constant, so that in the 
absence of a convection electric field, the Poynting flux vanishes. However, 
if we have a region between the neutral sheet and the end of the domain in 
which there are no particles, that is, where the source term for (28) 
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vanishes, the value of ■ da/dz will be constant, so It Is necessary to 
calculate the solution only where n(z) Is assured of being nonzero. 

Let the domain be divided Into two regions, 0<z<LandL<z<L + 0, 
and let n be zero In the latter region. We therefore need to Integrate (28) 
only In the region 0 < z < L, with the boundary condition at z ■ L, 
a + D5a/5z ■ 0. This mixed boundary condition Is easily Incorporated Into the 
finite difference, tri-diagonal algorithm (Roache, 1976). In this case, the 
total magnetic energy Is given by 

En. ■ 2^ i:/i< dz + 0 o 2(L) + (L . D)o^] (30 

Reflecting boundary conditions must be used in the particles at z » 0. 
This means that a particle which finds Itself at z < 0, the replacement z -z 
and Vj ^ -v^ is made. The other phase space variables are unchanged. These 
boundary conditions follow from the assumption of reflection symmetry. When a 
particle from the north passes through the symmetry plane, it Is replaced by a 
particle from the south side In which the z component of velocity Is 
reversed. This coupled with the fact that Q • 0, makes the z » 0 boundary 

A 

absolutely enerv^y conserving. Particles which exceed z ■ L are simply 
reinserted inside the domain, and the velocities are reinitialized to the E x 
B drift velocity. Although this particle boundary condition Is not energy 

conserving, very few particles cross the z L boundary. 


3, Results of Calculations, 

3.1 Two-Dimensional Model Results. 

The two-dimensional calculations seek to examine the question of how the 
earth's magnet otall and neutral sheet might be Inflated and respond to 
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changing external conditions, such as the Imposed convection electric field. 
The calculations also seek to explore whether turbulence would develop, and 
whether it would lead to the formation of X- and 0-type neutral lines. 

As a working hypothesis It Is assumed that a convection electric field Is 
present In a region accessed by cold particles of Ionospheric or magnetosheath 
origin. Panel a of Figure 1 shows the initial magnetic field configuration. 
This was generated as a contour plot of the vector potential. The 
magnetospheric z-axis Is upward, and the x-axIs, earthward direction. 
Increases to the right. In the region where the magnetic field Is uniform the 
X- and z-components of S, as It appears In (13a, b) Is set ot 0.5. The source 
of the x-component of the field Is a current sheet, which Is held constant in 
time, but has small random spatial perturbation. 

Contours of constant current density are exhibited In panel a of Figure 
2. It was necessary to introduce this current because the Interatlon 
procedure outlined in (2) requires a starting vector potential consistent with 
the initial particle distribution. This current sheet was assumed In order to 
Int.'oduce a kink in the magnetic field and the perturbation In the current 
sheet to excite any plasma Instabilities that might develop. It will become 
apparent that this ad-hoc current will be a small fraction of the total 
particle current. Some field-line curvature is necessary, in order that the 
E X B drift not be everywhere Inertial. It Is only in the region where the 
E X B drift is non- Inertial that ions will move relative to electrons. 

The planes shown in Figure 1 are defined by 129 x 129 grid points, which 
are loaded with 4 x lo^ particles in a band 129 x 89 grid points centered on 
the neutral sheet. The average density is 4 x 10^(128 x 88) » 3.55 particles 
per cell. This is the initial value of n in (16b). The initial particle 
velocities are aVl set to the local E x B drift velocity, such that they 


Initially contribute nothing to the current. The charge coupling parameter, 
a, In (15) Is set to a « 0.0655, which results In the asymptotic AlfvAn speed 
of V, - 3.040. The convection electric field Is 7.07, so that the asymptotic 
X and v-ccmponents of the E x B drift velocity are Ug ■ (7.07, 0, t7.07), with 
the z-components converging on the neutral plane. The Initial E x B drift 
energy of the plasma Is 2.179 x lo®, and the Initial energy In the magnetic 
field Is 6.U2 X 10^. Hence, in this model the E x B drift energy far exceeds 
the magnetic energy, for the boundary condition that the vector potential Is 
clamped on the z » 0 and z » L boundaries. The time-step for this model was 
At ® 0.025, and eight iterations were used each time-step to calculate Ay. 

Figure lb shows the magnetic field lines a short while later. The 
significant change is that the kink in the field lines Is sharper. Figure 2b 
shows the current density contours are also much narrower. The fixed current 
shown in Figure 2a Is unchanged, but the largest contribution to the current 
comes from the particles drifting Into the neutral plane. The current sheet 
appears thinner because the number of contours between maximum and minimum is 

i , 

constant. Independent of the difference between these values. 

Figures Ic and 2c show the field lines and current contours somewhat 
later. The current has bifurcated, and this Is reflected in the development 
of shock structures on either end of the neutral plane, which develop because 
the flow toward the neutral plane exceeds the phase velocity of a magnetosonic 
wave. This super-magnetosonic flow Is by virtue of the assumed cold particle 
distribution and very small magnetic field Intensity. It viould not be seen In 
an established plasma sheet and magnetotail in which particle thermal energies 
and magnetic field intensitites are much higher. The thickness of the shock 
is the order of resolution of the grid mesh. At the time corresponding 
Figures Id and 2d, all of the particles have been swept into the shock region. 
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and the external field begins to Increase. The asymptotic field lines appear | 

'f. 

somewhat curved In Figure Id because not quite enough Iterations have been 

taken to communicate the Information from the edge of the current sheet to the :| 

boundary of the domain. The feature of main Interest In Figure 2d Is the 

I 

development of multiple current sheets and the expansion of the current I 

carrying region. | 

Figure le shows the continued Increase of the external field and the ! 

continued thickening of the current carrying region. A satellite traversing 
the plasma sheet would observe multiple reversals In the x-component of the 
field. Figure 2e shows the development of current filaments and x-dependent 
structure. Plots of the z-component of the magnetic field at the original 

neutral plane show considerable short wavelength fluctuation. 

Finally, Figures If and 2f, at the end of the run show the development of 

turbulence In the plasma sheet region with some closed magnetic Islands. The 
run could not be carried on much longer because the current sheets would have 
expanded to the edge of the domain. The external x-component of the field has 
Increased to 3. .25, a factor of 6.5 over the x-component at the beginning of 
the run. This clearly Indicates that the action of a convection electric 
field on a cold particle population Is capable of significantly Increasing the 
lobe field. 

Figure 3 shows plots of the magnetic, kinetic and total energies during 
the run. The total energy began to increase at a more rapid rate only after 
all of the cold particles were swept Into the current sheet and the external 
field began to increase. Since the convection electric field was held 
constant, the value of Poynting flux at the boundaries would Increase In 
proportion to the boundary field. Also, It appears that the width of the 
current channel has a nearly steady Increase. 
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An important question is the fate of the turbulence seen in Figures If 
and 2f. The fact that the current sheet has successively broken into many 
smaller current sheets, and the fact that the irregularities have the same 
size scale on the x- and z-directions suggests that the- turbulent cells will 
not coalesce to form large scale structures of the type that would be 
identified with the development of a single or a few X-type neutral lines in 
the magnetosphere, or the development of the topology that emerged from the 
fluid calculations of Birn (1980). This conclusion is based on another run 
similar to that shown in Figures 1 and 2, except that the convection electric 
field was only half the value. This made it possible to carry the run on for 
a longer period of time. Turbulence also developed at about t » 15 into the 
run, and it was tracked until t ■ 20 at the end of the run. During the period 
the turbulence on the original center of the current sheet was observed to go 
through a cycle of growth, decay and regrowth, suggesting that the level of 
turbulence is stable. In this latter run, the diameter of the turbulent cells 
was the same as the width of the current sheets; however, the scale size and 
amplitudes were somewhat smaller than the scale size and amplitudes of the 
runs shown in Figures 1-3, with the stronger convection field. These results 
seem to be consistent with the results of fluid model simulations of Matthaeus 
and Montgomery (1981) of a sheet pinch which exhibited the development of 
small scale MHD turbulence. The turbulence did not lead to any enhanced 
energy dissipation. 

The most serious limitation of these runs is that the time span was not 
much longer than the characteristic ion gyroperiod. The possible consequences 
of this limitation will be discussed after the one-dimensional results are 
presented. Another point that should be kept in mind is that the magnetic 
field is very weak, so that the flows are super-magnetosonic. 


The model does seem to indicate a process for the inflation of the 

earth's magnetotail by a convection electric field. The model also suggests 
that, aside from low-level turbulence, a one-dimensional, z-dependent model of 
the neutral /plasma sheet may provide an adequate description. 

3.2. One-Dimensional Model 

An advantage of the one-dimensional model is that it is economically 

feasible to follow a system for a very large number of time-steps. Moreover, 
since the fields are calculated implicitly, no iteration is necessary, and 
action-at-a-di stance is possible. This makes it easy to follow rapidly 

propagating Alfv^n waves. 

The total model domain will extend from z ■ 0 to z » 256, but the 

particles will be confined to 0 < z •; L » 128. As explained in Section* 2.3, 

the z * 0 boundary is assumed to be the center axis of the current sheet with 

the boundary condition Q » 0. The vector potential is held fixed at z » 

A 

256. The model uses 8192 particles, divided initially into two populations. 
There is a hot population consistency of 2048 particles loaded into the region 

0 < 32. The particles are randomly distributed in z; their speeds are 

uniformly, randomly, distributed in the range 0 < v < 10/2; and the velocity 
directions are uniformly, randomly, distributed on the surface of a sphere. 
This population represents the plasma sheet. The charge coupling constant, oj 
of (28), was set to 3.9 x 10“^ in order that this population would support an 
asymptotic magnetic field of ■ 5. The other population consists of 6144, 
zero velocity particles which are spatially distributed according to 
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' ■ zng/32 z < 32 

■ Hg 32 < z < 112 

- 11^(127 - z)/15 112 < z < 127 

- 0 z > 127 

where Hq » 59.4. The z-component of the magnetic field Is ■ 1, so with 
n ■ 5.1, and a particle density of n^, the A1fv4n velocity Is » 10.6. Two 
separate runs were made with the one-dimensional model. One assuming a \ 

convection electric field, e , that Increases linearly from zero to 10 at 
t ■ 2w and thereafter remains constant. The other Is a nearly energy- 
conserving run with » 0. 

The results of the run with the convection electric field are displayed 
In Figures 4, 5, and 6, which show the particle number density, particle 
energy density and x-component of the magnetic field respectively at various 

stages of the run. The magnetic field lines are displayed In the contour 

plots of a(z) + OjjX In Figure 7, while Figure 8 shows a time history of the 
magnetic, kinetic and total energy. Figure 9 shows a plot of the final 
particle energy spectrum. 

The initial x-component of the magnetic field was computed self- 
censlstently with the plasma sheet particles, so If there were no z-component 
of the magnetic field and no electric field, the field and particle population 
would be in equilibrium. The stresses due to the addition of a z-component of 

I the magnetic field will cause an acceleration of the plasma in the x- 

I direction, while the y-directed electric field will cause a particle drift 

I 

i| toward the neutral plane. It can be seen that the cold portion contributes 

i 

I very little to the current. At t = 6 and 12, the effect of the convection 

electric field compressing the particle near the neutral plane can be seen. 
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This would be viewed as a thinning of the plasma sheet, a signature of 
substonm onset. The waves In the magnetic field profile are due to Alfv^n 
waves propagating outward from the plasma sheet, which are generated by 
Initial transients due to the fact that the Initial particle and field 
distributions were not chosen completely self-consistently. During this 
growth phase Figure 8 Indicates a slow Increase In magnetic energy and a much 
more significant Increase In kinetic energy. 

The middle portion of the run, from t * 30 to t « 70, saw an expansion of 
the plaimiR sheet, and a leveling off of the magnetic energy Increase, but a 
continuing and accelerating increase In plasma energy. Figures 4 and 5 show a 
merging of the cold and plasma sheet populations at about t » 50. In contrast 
to the plasma sheet expansion; Figures 6 and 7 show a continuing decrease In 
the thickness of the magnetic neutral sheet during this middle period. Also, 
the magnetic field profiles Indicate the existence of standing, shock-like 
features near the neutral sheet. 

The final phase of the run from t » 70 to 110, shows a continuing 
expansion of the plasma sheet, but now also shows a substantial thickening of 
the neutral sheet. At this time, all of the Ionospheric plasma has been swept 
into the neutral sheet and heated. The most significant feature of this later 
phase is the continuing increase in plasma energy and the comparatively rapid 
decay of magnetic energy. The magnetic field energy decays despite the 
continuing presence of the convection electric field. This indicates that 
maintenance of the tail field requires a continuing supply of cold particles 
to be swept into the neutral sheet region. 

Finally, Figure 9 shows the particle energy spectrum at the final stage 
of the run. Also shown for purpose of comparison is the energy spectrum 
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f(E)dE - 4.40 E/2e“^/® dE 

where e » 164. corresponding to a three-dimensional Maxwellian distribution of 
velocities that has the same average energy as the simulation energy 
spectrum. The most significant features is that the simulation spectrun is 

much more peaked than a thermal distribution. The spectrum is also 

character! zed by a long tail, containing relatively few particles that decay 
quite slowly, with the highest energy particles at 1819. 

Figure 10 shows magnetic, kinetic and total energy for a run which had 
the same initial phase space distribution of particles, but in which there was 
no convection electric field present. The result is that magnetic energy 
immediately begins to decay with the energy going into kinetic energy of the 
particles. This is an energy-conserving run, so the total energy should 
remain constant. The degree to which the total energy remains constant is a 

measure of the accuracy of the code. During the run, there was a 3.7% 
increase in total energy. The plots of density, energy density, and magnetic 
field during the early phase of the run showed little noteworthy differences 
between the run and the run with the convection electric field present, except 
that toward t = 40, the run with the convection electric field present began 
to show a thicker plasma sheet. Comparison of the two runs suggests that the 
presence of a convection electric field is necessary to keep the tail field 
from collapsing. During the field collapse the solenoidal electric field was 
the same order as the applied convection electric field in the previous 
example, =« 10. 

There are significant differences between the one and two-dimensional 
runs. The two-dimensional run exhibited strong shock formation. This is 
attributed to the fact that in the two-dimensional run the flow was super- 
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magnetosonic. The other major difference was the absorption of particle 
energy in the one-dimensional simulation. This is due to the fact that the 
one-dimensional simulation contained an energetic thermal population and that 
the one-dimensional run had a duration of several gyroperiods such that the 
cold population became heated and thermal ized. Another difference is that in 
the one-dimensional simulation, the E x B drift energy Increases by a factor 

of 26 between the asymptotic regions and the neutral sheet, while in the two- 

* 

dimensional simulation, there is only a factor of two increase in E x B drift 
energy. Moreover, in the two-dimensional run when the particles encounter 

shocks In which there Is a significant increase in the magnetic field, the E x 
B drift energy shows a decrease. 

4. Discussion and Conclusion 

The purpose of this paper has been to use one- and two-dimensional 

particle codes to model the formation and dynamics of the earth's 
magnetotail. The two-dimensional model was used to explore a process for the 
inflation of the earth's magnetic field Into a configuration like the 
magnetotail. The process assumed the presence of cold particles and a 
convection electric field and an initial kink in the magnetic field caused by 
the assumed presence of a sheet current. The simulation showed that the cold 
particles drifting into the region of the current sheet would expand the 
current sheet and eventually increase the magnitude of the external field. 

A question of interest is whether an arbitrarily small perturbation in an 
otherwise uniform field would give rise to the inflation exhibited in 
Figure 1. The firehose instability (Hasegawa, 1975) provides useful 
insight. The convection electric-field may be transformed away by going to an 

inertial frame moving at u = cEy/B^. The result is particle streaming 
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parallel to the magnetic field at Vj » cEyBj^/CBB^). Particles originating 
from opposite sides of the current sheet will stream In opposite directions, 
resulting In Interpenetrating Ion beams. If one assigns a parallel kinetic 
temperature to this distribution of 9^ , then the criterion for the 

firehose Instability says the kinks will grow when 

Bun 9 /B^ > 2 

" • 
tan^t|> > 2 (32 

where <|» is the angle of the asymptotic field lines make with the z-ax1s. For 
t a 45°, E » 1 mV/m and n » 1 cm"^, (32) is equivalent to B < 5.7y. 

Equation (32) can also be written in terms of the change in Bj^ across the 
current sheet, and it puts a lower limit on the size of a perturbation, given 
n, and Ey. Equation (32), when t is expressed in terms of Bj^ shows that 

there is a minimum perturbation, below which perturbations will not likely 
grow. These runs however indicate that a magnetotail configuration like those 
presently observed could evolve from a very weak field through a series of 
repeated injections of cold particles to build up a plasma sheet population 
which supports its extended magnetotail lobes. 

The other question addressed by the two-dimensional model is whether 
instabilities in the current sheet would likely develop which would lead to 
the formation of X- and 0- type neutral lines. The indications from the model 
runs is that fine-scale turbulence does develop and in one model it does lead 
to the formation of closed flux lines, but it does not seem to lead to the 
large scale x-dependence that would be identified with the formation of X- and 
0-type neutral lines in the magnetosphere. These conclusions must remain 
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tentative until more extensive runs can be made. The code was run under 
conditions of super-magnetosonic convective flow, which might exist in only 
1 ill 1 ted regions of the plasma sheet. It should also be run for longer time 
periods. Such runs must await extensiva revision to the program to Increase 
Its efficiency. 

The program also assumed a constant magnetic flux through the upper and 
lower boundaries of the simulation domain. This does not take account of the 
fact that as the magnetosphere Is Inflated, more and more field lines are 

pulled out Into the tall. This could change the magnetic flux through the top 
and bottom boundaries In a localized region represented by the simulation 
domain. A sufficient decrease in the flux might permit the formation of 

large-scale magnetic islands. 

The one-dimensional model also assimed a constant z-component of the 

magnetic field, and also x- independence of all quantities. In other respects 
it was more realistic, in that the initial magnetic field was supported by an 
initial distribution of hot particles, and that the convective flow velocities 
were less than the Alfv6n speed. Moreover, the system was followed for a 
large number of ion gyroperiods. A substorm was modelled by applying an 
external convection electric field and assuming the existence of a population 
of cold particles external to the plasma sheet population, representing 

particles of ionospheric origin in the lobes of the magnetotail. 

The plots of the energy density of the run shown in Figure 5, reproduce 
the well -documented thinning and re-expansion of *the plasma sheet. This 
thinning is attributed to the inward E x B drift of plasma sheet particles. 
This occurs whether or not a convection electric field is applied, because in 
the absence of a convection electric field, an induction electric field 
develops due to the collapse of the magnetic field. This induction electric 
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field also causes an Inward E x B drift. The subsequent expansion Is a result 
of more energetic particles streaming outward from the neutral plane. If the 
field-aligned component of the particle velocity Is large enough, the z- 
component of Its motion can exceed the Inward component of the E B drift, 
resulting In plasma sheet expansion. 

There also appears to be a thinning of the magnetic neutral sheet 

I 

followed by a subsequent thickening, which lags considerably behind the 
thinning and thickening of the plasma sheet. However, an observer stationed 
at one position with respect to the neutral plane would observe no clear trend 
In the variation of the magnetic field. The thinning would only show up as a 
more rapid reversal of the x-component of the magnetic field on traversing the 
neutral sheet. The thickening of the magnetic neutral sheet occurred after 
all of the cold particles had been swept Into the plasma sheet and 
thermal 1 zed. The thickening process also accompanied a decrease In magnetic 
energy. 

The decrease in magnetic energy at this time, in spite of the continuing 
presence of the convection electric field, suggests that a necessary condition 
for the maintenance of the tall "^ield is the continued influx of cold 
particles from the lobes of the magnetotail. This in turn implies the 

necessity of replenishing the lobe field lines with particles of magnetosheath 
origin. Figure 10, which shows the variation of magnetic and particle energy 
for a run in which there was no convection electric field, shows the immediate 
onset of magnetic energy decrease. This suggests that the presence of a 
convection electric field is also a necessary condition for the maintenance of 
the tail field. The fact that when both cold particles and the convection 
electric field are present, the magnetic field energy remains constant or 
increases, suggests that a necessary and sufficient condition for the 
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existence of the magnetotall Is the combined presence of the convection 
electric field and cold particle populations In the magnetotall lobes. 

The results of the simulation allow us to tentatively Identify the 
manifestations of a substorm In the magnetotall as simply a consequence of an 
Increase in the convection electric field. It appears that a quasl^steady- 
state field and particle configuration can be maintained with a modest 
convection electric field and a continuous supply of ionospheric or 
magnetosheath particles to the lobes. A temporary Increase In the convection 
electric field would bring about the thinning and re-expansion of the plasma 
sheet. There would be enhanced dissipation which would heat the plasma sheet 
population, but since there would also be other low energy particles in the 
process of being heated a satellite observer would not necessarily observe an 
Increase in plasma temperature. The lobe magnetic field could increase or 
decrease depending on the availability of low energy particles. The substorm 
decrease in magnetic energy observed by Fairfield at al. (1981) could follow a 
consequence of the depletion of the cold lobe particles. These processes take 
place without the occurrence of any large-scale plasma instability. 

The simulations suggested that the plasma sheet would continue to thicken 
whether or not a convection electric field is present, so there would be no 
way for the plasma sheet to return to its original configuration. Therefore, 
the complete cycle requires some process for the removal of plasma sheet 
particles. In the magnetosphere this can be accomplished through the 
convection of plasma sheet particles into the inner portions of the 
magnetosphere to become part of the ring current belt. 

Another important result is the fact that convection electric field or 
magnetic field energy can be readily converted into particle energy; and as 
predicted by Roederer (1977), this can occur in the absence of any magnetic 
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regiorv associated with a neutral line. Fluid codes, on the vither hand, 
require the special introduction of viscosity or resistivity in order to 
effect such a conversion (Birn, 1980; Matthaeus and Montgomery, 1980). The 
ready conversion of energy in ordered plasma flow to thermal energy is 
possible because as Swift (1977) and Wagner (1979) have shown, ions Interact 
non-adiabatically with the magnetic neutral sheet region where the ordered 
flow energy is converted into gyrational motion. This has the same effect as 
dissipation. 

The simulation time scale can be scaled to the magnetospheric time scale 
by the particle gyrofrequency. In the one-dimensional model, we have assumed 
a five-to-one ratio between the lobe and neutral sheet field strengths, which 
is not unreasonable. If we assume a lobe field of ISy, then a simulation time 
duration of 110 corresponds to a time of 383 sec if ions are assumed and 
6130 sec for O"^. The time assuming populations is shorter than actual 
substorm times. The time duration of the simulation could easily be increased 
by populating a greater portion of the magnetotail lobes with cold particles. 
One effect of a longer run would likely be a broader spectrum of particle 
energies than indicated in Figure 9. 

Distances can be scaled on particle gyroradius. A 1 keV ion is a 15y 
field has a gyroradius of 300 km. A numerical simulation velocity is 10, 
giving a simulation gyroradius of 2 grid units. The intial half-thickness of 
the simulation plasma sheet is 32 grid units, which would correspond to 
4800 km, somewhat less than observed plasma sheet thicknesses. An electric 
field of 10 units gives an E x B drift velocity of 10 units in the neutral 
plane, which corresponds to a velocity of 436 km/sec (1 keV) in magnetospheric 
units. In a 3y field, this corresponds to a convection electric field of 1.3 
mV/m. The equivalent number density of cold plasma can be obtained by 
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coniparlson of Alfvin velocities. In a magnetic field of 5 gyrofrequency units 
and a density of 59.4, the AlMn velocity 1s » 10.38, which corresponds to 
a magnetospheric Alfvin velocity of 458 km/sec. A number density of 0.53 cm“^ 
gives this Alfvin velocity In a 15 y field. Tbis number density Is well within 
the range reported by Anderson and Eastman (private communication). 

Finally, we wish to .assess possible effects of electrostatic fields 
arising from charge separation. Electrons and Ions will follow different 
trajectories In regions where the magnetic field Is non-uniform. As mentioned 
previously, we anticipatethat highly mobile electrons will be effective In 
maintaining charge neutrality along magnetic field lines. However, In regions 
where the magnetic field curvature Is large, Ions will stream perpendicular to 
field lines through strongly magnetized electrons. This could excite the 
modified two-stream Instability (McBride et al., 1972). This Instability can 
give rise to substantial electrostatic potential which would substantially 
change Ion trajectories. This is a subject currently under Investigation with 
a purely electrostatic code Involving both electron and ion dynamics. This 
process offers the intriguing prospect- of generating large potential 
differences across field lines. Such potentials may play a role 1n 
accelerating auroral electrons. 

Another possible Interaction Involving electrostatic fields Is the lower 
hybrid drift wave. This mode is excited by currents in the y-directlon, but 
plasma sheet density gradients act to destabilize the mode so that excitation 
occurs at lower current thresholds. Huba et al. (1977, 1978) point out that 
such Instabilities can give rise to anomalous resistivity and enhanced 
dissipation, perhaps sufficient to allow ion tearing mode Instability. 
However, the calculations presented here and individual particle calculations 
done elsewhere (Swift, 1977, Wagner et al., 1979) have shown that the cross- 
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tail electric field can result in the transfer of energy to gyrational motion 
of ions. This can have the same effect as dissipative term in the fluid 
equations, namely transfer of energy from ordered motion to internal degrees 
of freedom. Therefore, even if the model calculations reported here permitted 
excitation of the lower hybrid drift instability, it is- not certain that this 
would significantly effect the calculations. Moreover, Schindler and Birn 
(1978) have pointed out that when > 0 everywhere, and the tail field is in 
equilibrium with pressure gradients that it is stable against all modes. 
Although the conditions on the model do not involve equilibrium, the 
Schindler-Birn theorem suggests that the models used here might be stable 
against the tearing mode. 
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Appendix 


Details of Numerical Procedures to Solve Equation (21). 

The method of solution is similar to that described by Dycek and Dawson 
(1978). Equation (21) is Fourier transformed in x so that 

A(x, z) - I A(m, z) 
m 

where A(m, z) satisfies 

(■^ - sin^ -p - k^) A(m, z) » -p(m, z) (A1 

where N+1 is the number of grid points on [0 < x < L] and ax is the distance 
between grid points. 

We wish to solve (Al) subject to the boundary condition that the magnetic 
field must be asymptotically uniform. In order to impose this condition, we 
write (Decyk and Dawson, 1979) 

A(m, z) » a(m, z) + ?i(m, z) (A2 

where X(m, z) is the periodic solution, calculated from 

s:(m. z) - I ?(«., n) 

n 

where ?i(m, n) is given by 


n) - p(ni. n) f 4 , 4 sin^UsM * 

I AX AZ 


(A3 


where p(m, n) Is the double Fourier transform of p(r) in (21), and lz is the 
distance between grid points along the z-axis. 

The other term in (A2) satisfies the homogeneous part of (Al). Using the 
finite difference approximations to the second partial derivative, it can be 
shown that 


a(j, m) 


A„r*^ + 
m 


m 


-j 


(A4 


where j is the index denoting the grid point number along the z-axis, and 
, where n is given by 

2 

n ■ 1 + 2 sin^(iim/N) + k^Az^ (A5 

Az- 

For m > 1, the coefficients A^ and are determined by the condition that 
A(m, j) where j » z/Az, join smoothly onto solutions which vary as rgJ at j « 
0 and rg"J at j » N, where rg » + where is calculated as in 

(A5), but with set equal to zero. For the m » 0 solution to (Al), we 
impose the boundary conditions that 


and 


a'(0. 0) » -a'(0, N) 

I (A(0, 0) + A(0, N) - Ag 


(A6a 

(A6b 


where the primes denote the finite difference approximation to the derivative 
with respect to z. Equation (A6a) requires that the x-componentd of the 
magnetic field at z » 0 and L be equal and of opposite sign, while (A6b) 
requires that the value of AA/Az averaged over the boundary vanish. In the 


absence of the convection electric field, (A6b) specifies that there be no net 
Poynting flux across the boundary. 
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Figure 3* 



A plot showing the magnetic, kinetic and total energy during th< 
run shown in Figures 1. and 2. The energy Is in units of 0«5v 
where v is the particle velocity In simulation units. 
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Figure 4< Number denslCy profiles as a function of the distance from the 
symmetry plane, expressed In grid units. The density is 
expressed In particles per grid unit. The tic marks show 
densities In units of 100. 
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Particle energy density profiles for the run shown in Figure 
4. Conparlsons between the two figures make It possible to 
distinguish between regions occupied by cold and energetic 
particles* 


Figure 5 < 
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Figure 6. 


Profiles of the x-coiqponent of the magnetic field for the run 
depicted in Figures 4 and 5. The tic marks show a field 
strength of 5 in simulation units* 







figure 7. Contour plots of the vector potential A - a(z) + showing 
magnetic field lines for some of the profiles shown in Figure 


6. The bottom edge is at the neutral plane, while the top edge 


is at 2 - 128. 
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Figure 9. A plot showing the kinetic energy spectrum^ at the end of the 
run* The energy is in simulation units » V 2 V“* The smooth curve 
is the kinetic energy spectrum for a three-dimensional 
Maxwellian with the same number density and average energy* 
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Figure 10. A plot showing the time variation of kinetic, magnetic and total 
energy during a run similar to that In Figure 3, except that the 
convection electric field was assumed to be zero. 


